Matrix Exponential-Based Closures for the Turbulent Subgrid-Scale Stress Tensor 
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Two approaches for closing the turbulence subgrid-scale stress tensor in terms of matrix exponen- 
tials are introduced and compared. The first approach is based on a formal solution of the stress 
transport equation in which the production terms can be integrated exactly in terms of matrix 
exponentials. This formal solution of the subgrid-scale stress transport equation is shown to be 
useful to explore special cases, such as the response to constant velocity gradient, but neglecting 
pressure-strain correlations and diffusion effects. The second approach is based on an Eulerian- 
Lagrangian change of variables, combined with the assumption of isotropy for the conditionally 
averaged Lagrangian velocity gradient tensor and with the 'Recent Fluid Deformation' (RFD) ap- 
proximation. It is shown that both approaches lead to the same basic closure in which the stress 
tensor is expressed as the product of the matrix exponential of the resolved velocity gradient tensor 
multiplied by its transpose. Short-time expansions of the matrix exponentials are shown to provide 
an eddy-viscosity term and particular quadratic terms, and thus allow a reinterpretation of tradi- 
tional eddy-viscosity and nonlinear stress closures. The basic feasibility of the matrix-exponential 
closure is illustrated by implementing it successfully in Large Eddy Simulation of forced isotropic 
turbulence. The matrix-exponential closure employs the drastic approximation of entirely omitting 
the pressure-strain correlation and other 'nonlinear scrambling' terms. But unlike eddy-viscosity 
closures, the matrix exponential approach provides a simple and local closure that can be derived 
directly from the stress transport equation with the production term, and using physically motivated 
assumptions about Lagrangian decorrelation and upstream isotropy. 
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I. INTRODUCTION 

One of the most basic challenges in turbulence mod- 
eling is the need for closures for the fluxes associated 
with unresolved turbulent fluctuations. In the context 
of Large Eddy Simulation (LES), closures are required 
for the subgrid-scale (SGS) stress tensor [l|, 0]. Tradi- 
tional closures involve mostly algebraic expressions re- 
lating the stress tensor to powers of the velocity gradient 
tensor. More elaborate approaches using separate trans- 
port equations have sometimes also been employed, al- 
though these tend to be significantly more costly in the 
context of LES. Closures expressing the stress in terms 
of the matrix exponential function do not appear to have 
received much attention in the literature. The objective 
of the present work is to identify and discuss two sep- 
arate paths that lead to such closures. Both paths are 
based on the Lagrangian dynamics of turbulence, i.e. on 
an understanding of the evolution of turbulence as one 
follows fluid-particle paths in time. 

The use of Lagrangian concepts in turbulent flows has 
a long history [HJj] and, in recent years, has seen renewed 
interest for modeling 

H, H S H H • Among others, a new 
model for the pressure-Hessian tensor based on the re- 
cent Lagrangian evolution of fluid elements - the Recent 



Fluid Deformation (RFD) closure - has been proposed 
[ToL fill [l2j . In this approach, a change of variables is 
made expressing spatial gradients in terms of Lagrangian 
gradients (e.g. how does a variable at the present location 
vary if we change the initial position of the fluid particle 
at an earlier time). Then the assumption of isotropy is 
introduced for the Lagrangian gradient tensors. This as- 
sumption allows simpler isotropic forms to be used, and 
is argued to be justified based on Lagrangian decorre- 
lation effects. Deviations from isotropy at the present 
location for the Eulerian gradient tensors develop as a 
result of fluid material deformation along the Lagrangian 
trajectory. More traditionally, the Lagrangian time evo- 
lution of the stress tensor following fluid particles can 
be derived by taking appropriate moments of the Navier- 
Stokes equations. In this paper we examine both of these 
approaches to formulate models for the SGS stress tensor 
in turbulence in the context of LES. 

A description of small-scale structure of turbulence be- 
gins with the Navier-Stokes equations of an incompress- 
ible fluid of velocity u: 



du du . _. _ 
- = - + (u.V)u = -Vp- 



vV 2 u , 



(1) 



where d/dt stands for the Lagrangian material derivative, 
p the pressure divided by the density of the fluid and 
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v the kinematic viscosity. Because of incomprcssibility, 
the velocity gradient tensor Ay = dui/dxj must remain 
trace-free, i.e. An = 0, and the pressure field is the 
solution of the Poisson equation V 2 p = — AikAki- 

In the framework of LES, the SGS stress tensor is de- 
fined using the filtering approach [LH [Til l , 



Tij = UiUj - UiU 3 . (2) 

An overbar denotes spatial filtering at a scale A and is 
formally given by a convolution with a nonncgative, spa- 
tially well- localized filtering function Q (r) of characteris- 
tic size A, with unit integral J Q(r)dr = 1, namely 



u(x,t) = J g(r)u(x+r,i)dr 



The SGS tensor r enters in the dynamics of the filtered 
velocity u as it can be seen when applying the filtering 
procedure to the Navier-Stokes equations (Eq. |T])), 

Du du ._ __ _,_ _ 

— = + (u- V)u= -Vp + vV-u- V-t , (3) 

where D/Dt stands for the Lagrangian material deriva- 
tive with u as the advecting velocity, and p the filtered 
pressure divided by the density of the fluid. Because 
of incomprcssibility, the filtered velocity gradient tensor 
Aij = dUi/dxj must remain trace- free, i.e. An = 0, and 
the filtered pressure field is the solution of the respective 
Poisson equation V 2 J5 = — AjaAw — d 2 Tij / dxidxj . 

Next, we also consider the transport equation for the 
SGS stress tensor r fl3l . \\A [l5j , which follows from Eq. 
©: 

Dt 9t t 

= — + (u- V)t = -tA - Ar + $ , (4) 

Dt dt y 1 > w 

where the term $ = <& p +<I?„ V jT includes the pressure 
gradient-velocity correlation 



®p,ij = _ [uidjp - mdjp + Ujdip - Ujdip\ , 
the viscous term, 



^v,ij = f{Ui^ 2 Uj — UiV 2 Uj + UjW 2 Ui — UjW 2 Ui) 

and the generalized central third-order moment 



Ji 



ijk = UiUjU k - UjT ik - UiTjk - U k Tij - Ui UjU k ■ 

In fjn] it is shown that a formal solution for the stress 
transport equation may be obtained by integrating the 
production term exactly. This solution, suggested by 
[l9( but — to our knowledge — little pursued, will be shown 
to involve matrix exponentials. The developments pre- 
sented require some assumptions of Lagrangian isotropy 
and decorrelation, and some empirical evidence support- 
ing these assumptions is provided in §1111 based on re- 
sults from Direct Numerical Simulations (DNS). In §TV] 
the RFD closure for the SGS stress is developed. The 



resulting model is shown to be expressible compactly in 
terms of matrix exponentials as well. Differences and 
similarities between the RFD and transport equation so- 
lutions are discussed. In §Vl the matrix-exponential solu- 
tions are expanded for short times. The expansions allow 
to establish relationships to traditional eddy-viscosity 
and nonlinear closure models in turbulence. In §VII the 
matrix-exponential closure is implemented in a most sim- 
ple flow to illustrate its feasibility and cost. 



II. SOLUTION TO STRESS TRANSPORT 
EQUATION USING MATRIX EXPONENTIALS 

Equation[4]is of the form of the "time-dependent Lya- 
punov equation", if the tensor <&'s implicit dependencies 
upon the velocity fluctuations and the stress tensor were 
disregarded (in reality, $ depends upon small-scale ve- 
locity fluctuations and thus the full equation is highly 
non-linear and non-local). The formal solution of the 
Lyapunov equation in terms of matrix exponentials has 
been found useful in a number of other fields: principal 
oscillation pattern analysis |16j . mechanics of finite de- 
formations [l7j , and fluctuation-dissipation theorems for 
stochastic linear systems [HI- In the context of the SGS 
stress transport equation the solution at time t (start- 
ing from an initial condition at time t') may be written 
formally as follows 

i-(t) = H(t, t')T(t')H T (t, t')+ J H(i, s)*(s)H T (i, s)ds , 

(5) 

where 



DK pt^ = -A(f)H(t, t') and H(t', t') 



(6) 



To our knowledge, this approach to solve the stress equa- 
tion in RANS closures was first suggested in the turbu- 
lence literature by [l9[ (see equation (4.4) in this ref- 
erence). For the general case of time- varying veloc- 
ity gradient, we note that the auxiliary matrix H(t, t') 
can be written as a time-ordered exponential (see Refs. 
[3 HH for background on this basic matrix function) 



H(i, t') = Texp 



- / A(s)ds 



Equation ([5|) illustrates clearly the distinct roles played 
by the production term and the contribution given by 
Evaluation of Eq. [5] requires the knowledge of the time 
history of A(s) as well as accurate closures for 3>(s) along 
the fluid history t' < s < t. 

As a next step, one may consider the special case 
in which the velocity gradient is considered to be con- 
stant between the initial time t' and t, and set equal to 
(e.g.) A(t) and simply denoted by A. For this approx- 
imate situation, the solution of Eq.Q may be written 
as an ordinary matrix exponential H(f, t') = e - '* - * ) A , 
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where the matrix exponential is defined in the usual way 
e B = J2n°°o Bn / nl To simplify further, consider Eq. |5] 
for the case = 0, i.e. now retaining only the production 
term. This step eliminates the important isotropization 
effects of pressure-strain and also the nonlinear diffusion 
effects of the transport terms. While clearly missing im- 
portant physics, it is still instructive to observe that this 
simplification allows the solution (Eq. O to be written 



T(t) 



,-(t-t')A 



T(0 



-(t-t')A' 



(7) 



At this stage it is conceptually advantageous to make 
connection with Refs. [H|,[23], 24], where it is proposed to 
use conditional statistics to capture the relevant statis- 
tics of the SGS stress. For example, in Ref. [23) . it is 
shown that the least-square-error best estimate for the 
SGS stress is of the form of a multi-point conditional av- 
erage, namely (r^ | Ui, U2, ...,Ujv)- The multi-point con- 
ditioning variables {ui,U2, ...,uj\r} are, in principle, con- 
stituted by the entire (./V-point) resolved velocity field at 
scale A. To simplify the conditioning, one may limit the 
information to the past time-history of the local veloc- 
ity structure. In particular, a good choice that captures 
much of the local dynamics in a Galilean invariant fash- 
ion is the Lagrangian past history of the filtered velocity 
gradient tensor A. The dependence on the Lagrangian 
time history along a fluid particle advected by the filtered 
resolved velocity field is thus assumed to be described by 
A(s) with s < t (here and below, the dependence of A 
on spatial position is omitted for clarity). According to 
these ideas, we define a quasi- optimal SGS stress ten- 
sor T^°\t) as the conditional average T^°'(t) = (t(£)|A). 
Noticing that the matrix exponential prefactors entering 
in Eq. [7] arc some deterministic functions of the velocity 
gradient tensor itself, they thus can be taken out from 
this conditional average and we get the following stress 
tensor: 



(*) 



,-(t-t')A 



[r(t')\A) 



,-(*-t')A' 



(8) 



With this expression, the closure problem has been 
changed from requiring a model for the local stress ten- 
sor at time t to requiring a model for the conditional 
average of the 'upstream initial condition' at time t' < t. 
The initial condition needed is a symmetric tensor. In 
the absence of additional information, the simplest as- 
sumption is to postulate that this conditionally averaged 
'upstream' stress tensor is isotropic, namely 



,(*')|A 



1 



-<r fcfe (t')|A) 6ij. 



(9) 



The magnitude of the tensor is proportional to the 
trace of the SGS tensor and has units of squared ve- 
locity. The assumption of isotropy may be justified if 
r(i') and A(£) become more and more de-correlated as 
the elapsed time t — t 1 grows, then no locally strong and 
statistically preferred direction should exist. This step 
introduces a characteristic dccorrclation time-scale r a , 



and t — t' will be chosen to be of the order of such a 
decorrelation time-scale. Clearly, one must also assume 
local isotropy to hold for the statistics, and this is jus- 
tified from the usual arguments in turbulence when A 
is sufficiently small compared to the integral scale. In- 
cidentally, it is expected that a decorrelation between 
r(t') and A(i) may occur due to pressure effects, tur- 
bulent diffusion, etc. Some numerical evidence for such 
decorrelation and isotropization is provided in the next 
section. 

The trace of the conditional SGS tensor, (Tkk{t')\A(t)) , 
must still be specified. The simplest option that is con- 
sistent with a local evaluation of velocity and length- 
scales is to choose a factor proportional to A 2 |S| 2 , where 

S = (A + A T )/2 is the filtered strain rate tensor and 
|S| = (2SijSij) 1/2 . Finally, replacing into Eq. [8] with 
t-t' = r a , we get 



Jo) 



D A 2 |S|V 



i„A„-f„A 



(10) 



where the parameter c oxp is unknown and may be ob- 
tained by empirical knowledge, or by generalizing the 
dynamic model [25l |. 

For completeness and clarity, we remark that the ma- 
trix exponential solution may equivalently be obtained 
by solving the linearized equation for a turbulent fluc- 
tuation that only keeps the Rapid Distortion term from 
the large-scale velocity field, and neglects all other ef- 
fects. That is to say, we solve formally the equation 



DM 



,Aij~ using the matrix exponential function. 



The solution is then multiplied by its transpose to form 
?4(£)iij-(i) which is then averaged over the fluctuating ini- 
tial condition u^(i')itj-(i') (conditioned on a constant A). 
The averaging of the term it^(£')u^(i') yields the initial 
('upstream') stress tensor r(£'), and with the conditional 
averaging, an expression equivalent to Eq. [8] is obtained. 
This is similar to the equivalence between solving the 
equation for co- variances or for the fluctuations and then 
averaging, as noted in the context of stochastic linear 
systems in [18 1. 

POP represents a closure for the SGS stress 



Equation 

expressed in terms of matrix exponentials instead of the 
more commonly used algebraic closures [26j . In section 
HVl a connection is noted between the expression Eq. (fT0|) 
and a physical closure for the subgrid stress tensor based 
on the recent fluid deformation closure in the Lagrangian 
frame. 



III. EMPIRICAL EVIDENCE OF LAGRANGIAN 
DECORRELATION AND ISOTROPY FROM DNS 

In order to verify whether the decorrelation and 
isotropization of conditional averages of SGS stresses oc- 
cur in turbulence, we analyze a DNS dataset of forced 
isotropic turbulence. The simulation is conducted using 
a pseudo-spectral method in a [0, 27r) 3 box. 128 3 grid 
points arc used. Fourier modes in shells with |k| < 2 are 
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forced by a term added to the Navicr-Stokcs equations 
which provides constant energy injection rate €f = 0.1. 
The viscosity of the fluid is v = 0.0032. Data is collected 
after the simulation reaches statistical steady state. Note 
that r(t') is the SGS stress at a previous time t' and at 
the spatial location occupied by the fluid particle which 
is at the position x at time t (i.e. X(i';x, t), in the no- 
tations of jjlVp . According to the transport equation for 
r(t) (Eq. H}, the fluid particle is advected by the filtered 
velocity field. Thus, the position of the fluid particle at t' 
is found by backward particle tracking starting from end- 
time t in the filtered velocity field. To perform backward 
particle tracking, the filtered velocity and SGS stress 
fields are calculated and stored at every At = 0.009, 
corresponding to 1/20 of the Kolmogorov time scale. A 
Gaussian filter is used with filter scale A = where 
r\ is the Kolmogorov length scale. In order to quantify 
isotropy as function oit — t', the ratio of off-diagonal to 
on-diagonal tensor elements of the conditional averaged 
SGS stress at decreasing previous time t' is computed. 

According to the derivation, the averaging must be 
conditioned on a particular value of the resolved veloc- 
ity gradient A(t). There are a large number of pos- 
sibilities, since A(t) has 8 independent elements. As 
representative of an important class of velocity gradi- 
ent structure, we choose to consider regions where the 
A(t) is such that it has a large shear in one direction, 
whereas all other velocity gradient tensor elements are 
weak. We choose a particular shear direction, "12", and 
define Ei2(i) to be a "high-12-shear" events that oc- 
cur at time t. These events are defined here as those 
points where 1) A\2(£) > A rms , i.e. large and posi- 
tive 12-shear, 2) |Ay(t)| < ^4 rms for other off-diagonal 
components = (1,3) and = (2,3), and 3) 

|Ay (t)| < A ims /v2 for the diagonal elements i = j. The 
gradient rms ^4 rms is defined as A rms = (AijAij) 1 ' 2 . This 
definition allows a sufficiently large number of events to 
be counted and thus help in reaching statistical conver- 
gence. With this definition of a conditioning event, we 
calculate the isotropy factor I(t — t') according to: 



0.15 



I{t-H) = - 



(T 12 (f')|E 12 (f)) 
(l/3)<7* fc (t')|Ei2(t))' 



(11) 



I(t — t') monitors the isotropization of the SGS stress 
associated with large "12 shear events", i.e. a particular 
anisotropic condition in the large-scale velocity gradient 
tensor. Since the turbulence is statistically isotropic, sim- 
ilar results are expected if the other two shear component 
of Tij, namely T13 and T23, had been chosen instead of T12, 
under conditioning based on events Ei3(t) or E23(i), re- 
spectively. 

Backward particle tracking starts from spatial loca- 
tions where the conditions in Ei 2 (t) are verified, at time t. 
At each time t 1 < t, the particle locations are calculated 
from the stored filtered velocity fields using a second- 
order Adam-Bashforth scheme. The filtered velocity and 
SGS stresses at the particle locations are interpolated 
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FIG. 1: Decay of the anisotropy factor / (Eq. 1 1 1 j) as func- 
tion of normalized time lag (t — t')A lms measured in DNS of 
isotropic turbulence. 




FIG. 2: Decay of the correlation coefficient p (see Eq. I12[) 
between mit — t') and —A\2{t) as function of normalized time 
lag [t — t')A Iln s measured in DNS of isotropic turbulence. 



from the stored fields using 6th order Lagrangian inter- 
polation. The conditional averages are then found by 
averaging over all tracked particles. Statistical sampling 
is increased by averaging over the trajectories starting 
from several different end times t and also over the other 
two 13 and 23 off-diagonal elements (in both Tij(t') and 

Ey(i)). 

The resulting ratio I(t — t') is plotted in Fig. [T] as a 
function of the normalized time lag (t — t')j4 rms . It is 
evident that the conditional average of the SGS stresses 
becomes more isotropic as the time lag increases and I(t— 
t') crosses zero at about 0.7 eddy turn-over times, namely 

t ~ 0.7A xma . Then there is negative undershoot to about 
negative half of the initial value, before it is relaxed to 

around zero (the isotropic value) at about t ~ 6A rms . 
The undershoot below zero is an interesting trend and 
understanding the physics of this behavior would be an 
interesting goal for future studies. 
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As additional evidence for the Lagrangian time decor- 
relation between stress and large-scale velocity gradient, 
in Fig. [2] we show the correlation coefficient between 
T\ 2 (t — t') and —Ai 2 (t), namely 



(T 12 (t-f)A 12 (t)) 

(r 12 (t-t>r)(A 12 (ty 



(12) 



The correlation is near 15% at zero time-lag (similar 
to the correlation coefficient between SGS stresses and 
strain-rate tensor often quoted in a-priori studies), but 

then decays to nearly zero at times around t ~ 2A Il ^ la . 
Taken together, the DNS analyses thus provide evidence 
for the isotropy assumption on (Tij(t')\A(t)}, as long as 

t - t' > r a with r a = A ims , 

Note that due to the cost of storing the entire sim- 
ulation for backward particle tracking, only moderate 
Reynolds numbers were considered in the analysis. The 
forcing length scale has been estimated to be about 50 
times the Kolmogorov length scale, 77, and the viscous 
effects begin to significantly damp the motions at scales 
of about IO77 and smaller. Therefore, using A = 15r?, 
there may be some effects from the forcing and viscous 
scales on the results. However, the observed tendency 
towards isotropization is expected to become more, not 
less, prevalent at higher Reynolds numbers. We point out 
that opportunities for much more in-depth future anal- 
yses of such issues are provided by the availability of a 
turbulence database at higher Reynolds number [29| (al- 
though this database could not be used for the present 
data analysis due to the fact that it does not yet contain 
sufficiently efficient means of filtering the data). 



IV. STRESS TENSOR MODEL BASED ON THE 
RECENT FLUID DEFORMATION CLOSURE 

This alternative approach is based on relating the SGS 
stress tensor to small-scale velocity gradients. To beg in, 
one may recall the multiscale expansion [2?], HH, l30t | in 
which, among others, the exact subgrid stress (Eq. ^ is 
written in terms of u 5 , the velocity field coarse-grained at 
a scale 5, but still with 6 << A, i.e. containing significant 
contributions from sub-grid scales. One may then define 



the approximated stress tensor r 5 



and 



naturally r = lim^o T ■ 

Consistent with the Kolmogorov phenomenology, as 
argued formally in j3lj . and also as used in various a- 
priori analyses of experimental data (see e.g. [28|) the 
SGS stress is relatively local in scale, stating that the 
leading terms entering in its development are given by 
the coarse-grained velocity at the resolution scale A and 
including also the next range of length-scales between 
5 « A//3 and A (e.g. (3 ~ 2). As a consequence, one may 
use the approximation T%j m Tj*~ A . Furthermore as- 
suming that u 5=A /^ is sufficiently smooth over distances 
A (or using the 'coherent subregion approximation' [3l|). 



a Taylor expansion of and evaluation of the filtering 
operation at scale A in Eq. [5] leads to 



C 2 A 2 



dx k dx k ' 



(13) 



One observes that similarit y- ty pe models such as the 
standard nonlinear model (2l.l32l] correspond to using the 
gradient of the large-scale velocity field (6 = A or /3 = 1). 
Nevertheless, it is the case /3 > 1 which is physically rel- 
evant since the true SGS stress includes scales smaller 
than A. However, for (3 > 1, the expression [T2] does 
not constitute a closure since then u s contains sub-grid 
motions that are not known at the LES filter scale A. 

As in @] and Chevillard & Mencveau (2006 - CM06 
from here on), a Lagrangian label position X is em- 
ployed to encode the time-history information. Using the 
two-time formulation of Q, the label positions X(i'; x, t) 
satisfy dX/dt' = u(X(i'),*') with X(i) = x. Thus 
X(i';x, t) represents the position X at a prior time t' of 
the fluid particle which is at position x at time t. Making 
the Eulcrian-Lagrangian change of variables also used in 
CM06 leads to the following expression: 



C 2 A 



2 dX p dX q du\ 



di 



dx k dx k dX v dX a 



(14) 



All terms in this expression are strongly fluctuating vari- 
ables. But, as in prior section, the most relevant informa- 
tion is retained by the conditional averaged expression. 
We propose the same conditional averaging based on the 
time-history of the velocity gradient tensor along the past 
fluid particle trajectory. Therefore, combining the con- 
ditional averaging and the change of variables one may 
write 



' dX p dX q duf du s 3 



dxk dxk dX p dX q 



A(s);t' <s<tj, 
(15) 

where the dependence of stress (t) on current position 
x is understood and not indicated to simplify the nota- 
tion. The Jacobian matrix Gij(t',t) = dXi(t';x,t)/dxj 
satisfies (see for instance [i~7^ 

D t G(t', t) = -G(t',t)A(t) with G(t, t) = I , (16) 

where I is the identity matrix. Thus, 



G(t',t) = Texp" 



A(s) ds 



is expressed as an "anti-time-ordered exponential" , with 
matrices ordered from left to right for increasing times 
[20I I21I . 33] . The only difference with the matrix function 
H(t', t) of the preceding section (Eq. ([6])) is the sense of 
time-ordering. 

Since the deformation gradient tensor G p k = dX p /dxk 
is a deterministic function of the past velocity gradient 
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history, these tensors can be taken outside the conditional 
averages in Eq. 1151 So far one can thus write 



dx k dxi 



with Yi 



dU S ; Oil* , _ 



ijpq 



dX p dX q 



A(s);t' <s<t 



where Y is a 4th rank Lagrangian gradient tensor. At 
this stage, it is now possible again to invoke Lagrangian 
isotropy, following the approach of CM06 and of the 
preceding section. It is assumed that the tensor Y is 
isotropic due to loss of information caused by turbulent 
dispersion, past pressure effects, etc. if t — t' is long 
enough. Under the Lagrangian-isotropy closure assump- 
tion, one may write 



Yi 



ijpq 



^4 Sij 5pq 



B SipSjq 



C'SiqSjp. 



(17) 



While individual realizations of a small-scale gradient 
tensor in turbulence are of course not isotropic, statis- 
tical moments such as the conditional average can be 
more justifiably approximated as isotropic. The isotropy 
assumption states that the rate of change of turbulent ve- 
locities 6u s (x, t) (at the present location and time (x, t)), 
with respect to changes in past locations of the fluid par- 
ticles at time t', is insensitive to orientation of <5X. This 
appears to be a plausible postulate, if sufficient time has 
elapsed, i.e. if t — t' is sufficiently large for decorrela- 
tion to take place. In the preceding section we used data 
from DNS to test the accuracy of such a de-correlation 
and 'Lagrangian isotropy assumption' in a closely related 
context (directly based on the stresses rather than small- 
scale velocity gradient statistics). Still, it is important to 
recognize that this step is introduced here as an 'ad-hoc' 
closure assumption and no claim is made that this is a 
formal step with controlled errors. 

While the assumption of isotropy eliminates the de- 
pendence of Yijpg upon A, the latter still affects the Ja- 
cobian matrix Gij = dXi/dxj that enters in the closure 
for the SGS stress. Next, we focus attention only on 
the trace-free part of the modeled SGS stress tensor, i.e. 
we subtract the trace of the stress. And, noticing that 
the 'right' Cauchy-Green tensor dkX p dkX q is symmetric, 
only the unknown B' + C enters in the resulting 'quasi- 
optimal' model for the deviatoric part of the stress (su- 
perscript od) model T^- . Dimensionally, the parameter 
B' + C has units of inverse time-scale squared, and de- 
pends upon the turbulence statistics down to scales 5. 
For a fixed ratio A./S, and with both scales in the iner- 
tial range, for simplicity we assume that the parameter 
B' + C follows, as in the prior section, 'Smagorinsky 



scaling', i.e. B' + C 



One thus obtains 



Sod) 



(0 



1 DXrn dX n 

3 dx k dx k 



where the parameter c cxp — C2 ■ c, in a similar fashion as 
in the preceding section, is unknown and may be obtained 
by empirical knowledge, or by generalizing the dynamic 
model H|. 

As a final step, the Recent Fluid Deformation (RFD) 
approximation is used (CM06) in which the time- varying 
velocity gradient A(s) between t' and t is approximated 
with a constant value (e.g.) equal to its value at t 
and denoted by A. The initial condition for the fluid 
deformation (when the deformation gradient tensor is 
assumed to be the identity), is prescribed at the time 
t' < t. The solution to Eq. [TH] can then be written as 
G(t',t) = e~(*~* ■ )A . Note that in this approximation, 
H(£, t 1 ) = G(t' , t), since the sense of ordering of the ma- 
trix products is no longer significant. 

The next step is to replace the solution for G(t' , t) into 
Eq. [TU And, as was done in the preceding section, to 
assume that a characteristic de-correlation time-scale r a 
has elapsed between the time where the initial isotropy 
assumption is justifiable and the current time when the 
stress closure is required. This means replacing the initial 
time t' with t — r a . Finally, the closure for the deviatoric 
part of the stress reads 



Sod) 



--cxp 



A 2 |S| 



(19) 



where all quantities are evaluated at (x, t) . It is imme- 
diately apparent that this closure is equivalent to the 
formal solution developed in the previous chapter (see 
Eq. (HUD. 



V. EXPANSIONS 

In preceding sections it has been shown that a matrix- 
exponential closure for the deviatoric part of the SGS 
stress tensor may be written as in Eq. As a next 

step, the behavior of this closure is explored when r a is 
small enough so that the norm of r Q A is much smaller 
than unity. Then e"^ psI - t q A + (l/2)(r a A) 2 + .... 
Up to second order one then obtains 



,A 2 |S| 



2 t„ S 



aa t + 1(a 2 + (a t ) 2 



(20) 



(18) 



Crow, in Ref . [19( , derived essentially the same result (see 
his eq.(5.3)) but with unspecified coefficients obtained as 
moments of his memory kernel and an additional term 
proportional to the material derivative DtS. It is imme- 
diately apparent that if the time-scale r a is chosen as 
r a = |S| _1 , then the first term is the standard Smagorin- 
sky model with c oxp = c 2 . (where c s is the Smagorinsky 
coefficient). Furthermore, the second term, the term in 
the square parentheses, is of the form of the 'nonlinear 
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model' [3, [H, with a prefactor c oxp A 2 . Two differ- 
ences with the standard 'nonlinear model' are apparent, 



however. The first is that if c c 



s?, then as coefficient 



of the nonlinear term this is significantly smaller than the 
coefficient for this term normally mentioned in the liter- 
ature (which ranges typically between 1/12 to 1/3). The 
second difference is the presence of the additional term 

2 T \ 

A + (A ) 2 ) /2. To make connections with standard 

non-linear models used more often in RANS (e.g. [3~i|). 
the velocity gradient is decomposed into symmetric and 
antisymmetric parts, A = S + $7. The result is (again 

with Ta = isr 1 ) 



-2c 2 A 2 |S|S 



c 2 A 2 



s 2 



-fns-sn) 

2 v ' 



(21) 

It is interesting to note that the expansion including the 

2 T 

term (A + (A ) 2 )/2 cancels exactly the $7 $7 part that 

is included in the standard non-linear model A A . For 
detailed a-priori studies of the various decompositions of 
the velocity gradient and non- linear terms see [HI]. 




FIG. 3: Radial kinetic energy spectra of forced isotropic tur- 
bulence from LES using the matrix-exponential closure of Eq. 
[19] with c cxp = (0.1) 2 . Solid line: 7 = 1, dash-dotted line: 
7 = 2, and dashed line: 7 = 0.5. Dotted line: universal 
Kolmogorov spectrum E(k) = 1.6ej //3 fc~ 5/ '' i . 



VI. MATRIX EXPONENTIAL CLOSURE IN 
LES OF ISOTROPIC TURBULENCE 

The expansion introduced in the last section is for- 
mally valid only for small values of the norm of T a A. For 
more realistic larger values, the expansion may be inaccu- 
rate and many additional higher-order terms are needed. 
They can all be expressed in terms of expansions into in- 
tegrity bases (26[, but it is in general difficult to obtain 
the coefficients of the expansion. Instead, it is proposed 
here to utilize the matrix exponential directly in simu- 
lations. Since the exponential involves the full velocity 
gradient tensor, it appears more natural to choose the 
time-scale t q according to t q = j(Aij Aij) -1 / 2 = 7|A| _1 
instead of using | S | 1 - The parameter 7 is an empirical 
coefficient of order unity. 

As a first test, LES of forced isotropic turbulence is 
performed. This flow is the simplest possible test-case 
and it is used here simply to determine whether simula- 
tions using the matrix-exponential based closure are nu- 
merically stable yielding realistic energy spectra, and to 
ascertain the associated computational cost. The gener- 
alization to dynamic versions and tests in more complex 
flows will be left for future investigations. The simulation 
uses the same pseudo-spectral method as was used in the 
DNS outlined in gTTJ with same grid resolution, forcing 
scheme and time step size. Dealiasing is performed by 
zero-padding according to the two-third rule. The vis- 
cosity of the fluid is v = 0.000137. The subgrid-scale 
model implemented is given by Eq. [Inland c cxp = (0.1) 2 
is chosen (dynamic versions [25j] of this model to deter- 
mined c cxp can be developed in the future). To specify 
r a , the values 7=0.5, 1 and 2 are tested (a dynamic ap- 
proach of determining 7 could also be developed) . In the 
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FIG. 4: Longitudinal derivative skewness coefficient 5 as func- 
tion of simulation time. Lines are the same as in Figure [3] 
tl ~ 6 is the integral time scale. 



pseudo-spectral scheme, the modeled SGS stress is eval- 
uated in physical space and is made trace-free (this only 
affects the effective pressure, not the dynamics) before 
computing its divergence in Fourier space. 

The matrix exponentials are evaluated using truncated 
Taylor expansion with scaling and squaring [36l | . Specifi- 
cally, we need to evaluate exp(B), where B = — 7A/|A|. 
For a matrix C in general, the Kth order truncated 
Taylor expansion uses matrix polynomial Tr-(C) = 
S^=o C n /n! to approximate exp(C), incurring an error 
bounded by || C || A ' +1 /{[l- || C || f{K + 2)}{K + 1)!}. 
The error decreases with the norm of the matrix C. 
Therefore, to evaluate exp(B), we first define C = B/2- 7 , 
where the value of the integer j is chosen to ensure 
|| C ||< 1/2. exp(C) is then approximated by Tr-(C) 
and finally exp(B) is given by [Tx(C)] 2J . The cost of 
calculating Tr:(C) is reduced by using Caylcy-Hamilton 



8 



theorem to express C n (n > 2) in terms of I, C, C 2 , and 
the invariants of C. Choosing K = 7 , we obtain the 
following equation for T7(C) with an error smaller than 
lCT 8 : 



T r (C) = C I + CiC + C* 2 C 2 



(22) 



where 



QcRc Qh^c 



1 



4c_ 
3! 
Qc 
4! 



5! 
Rc 
4! 

^£ 
5! 



Qh 

5! 
6! 



7! ' 
2Q c i?c 
6! 

2QcRc 
7! 



7! 



Here Q c = -Tr(C 2 )/2 and i? c = -Tr(C 3 )/3 are the 
two non-zero invariants of C (note that Tr(C) = 0). 
In terms of cost, the above algorithm uses about (1 + 
j)N 3 + 5N 2 + 2N + 37 flops to calculate exp(B) when 
B is given, where N is the dimension of the matrix. In 
our tests j = 1 + floor (log 2 7), so j = 1 when 7=1 and 
the cost is estimated at about 140 flops for each stress 
evaluation. This can be compared with the single matrix 
multiplication needed for the nonlinear model, which is 
about N 3 ~ 30 flops. Overall with this closure, our code 
took about twice as long to run as compared to using the 
mixed model. 

Simulations were initialized with random Fourier 
modes and evolved until statistical steady state was ob- 
tained. No numerical instabilities were observed for the 
three parameter cases considered (c oxp = 0.01, 7 =0.5, 
1 and 2). In Figure [3] the energy spectra obtained from 
the three simulations as averaged in the time interval be- 
tween one and three large-eddy turnover times are shown. 
Figure [5] shows the time-evolution of the derivative skew- 
ness coefficient S = ((diTii) 3 ) / ((<9i?Zi) 2 ) 3 / 2 . As can be 
seen, the case 7=1 appears to yield physically meaning- 
ful results, which can be compared with the well-known 
results of the Smagorinsky model and the mixed model 
(see, for example, j4Q[). But, there is clear dependence on 
the parameter 7. The skewness coefficient quickly drops 
to values near —0.3 for 7 = 1 and —0.36 for 7 = 2. These 
are realistic values for filtered turbulence [33] . The skew- 
ness values for 7 = 0.5, on the other hand, appear to 
be too close to zero, consistent with some pile-up of the 
spectrum at high wave- numbers. 



VII. DISCUSSION 

A new closure based on matrix exponentials and as- 
sumptions about short-time Lagrangian dynamical evo- 
lution has been proposed. Matrix exponentials as for- 
mal solution of the stress transport equation provides 
interesting insights into the effects of the production 
(gradient-stretching) term. Historically, in the context 
of RANS modeling using additional transport equations, 
the (closed) production term has justifiably not been the 



focus of attention in the literature. In LES, however, due 
to practical constraints 'algebraic' closures are most often 
preferred. The present approach shows that the effects of 
production in the context of such closures may be taken 
into account directly based on an exact solution of the 
stress transport equation. A central step in the present 
approach is to use isotropy for the 'upstream' initial con- 
dition. Evidence for such isotropization of initial con- 
dition, given present large-scale velocity gradients, has 
been obtained using a DNS database. Implementation of 
the closure in LES of forced isotropic turbulence yielded 
good results. The computational cost is significant, but 
it is not prohibitive. Since our code with this model took 
about twice as long to run as with a traditional algebraic 
closure, LES with this model at a resolution of N 3 has 
similar CPU cost as LES with a traditional model run at 
a resolution of {2 1 ' A N) 3 - (1.27V) 3 . 

It is crucial to stress that the additional, more subtle 
physics of the remaining terms in Eqf4] (pressure effects, 
turbulent diffusion, dissipation, etc..) are, in general, 
unlikely to be well-represented by the simple assump- 
tion of 'upstream' isotropy. In addition, non-equilibrium 
conditions in which A varies quickly along the particle 
trajectory are not included in the closure as written in 
Eq. 1191 in which the velocity gradient is assumed to 
have remained constant over a time-scale r a . To explore 
non-equilibrium effects, the full time-ordered exponential 
function must be used, although this would still leave out 
the non-equilibrium effects of $ . To compare the present 
approach to other closures will require more in-depth 
testing in more demanding, complex flows (e.g. where ef- 
fects of anisotropy, non-equilibrium, and pressure-strain 
correlations are expected to be important). 

It is also instructive to consider the case of two- 
dimensional (2D) turbulence. Nothing in the closure 
strategics pursued here limits their application to space 
dimension three, at least nothing very obvious. How- 
ever, the expansions (Eqs. I20l2ip show that this is not 
likely to be a qualitatively good closure for space dimen- 
sion two, since one there expects an effective "negative 
eddy-viscosity" corresponding to inverse energy cascade 
[39[ ] . It is thus worth reflecting on some of the reasons for 
the inaccuracy of the closures in 2D, since this may help 
pinpoint potential shortcomings in 3D as well. First, it 
is known that the 2D inverse cascade is less local than 
the 3D forward cascade, with most of the flux coming 
from triadic interactions for a scale-ratio (3 = 4 ~ 8. 
[HI, [3^]. However, the starting point of the RED clo- 
sure, Eq. lTHI is not accurate for (3 ^S> 1. To get a qual- 
itatively reasonable alternative at (3 substantially larger 
than 1 — which involves only first-order gradients — one 
must instead use something like the "Coherent Subrc- 
gions Approximation" of [391 ]. On the other hand, the 
starting point of the closure approximation in Section [TT1 
the stress transport equation [H is exact in 2D just as in 
3D. The failure of the closure procedure in 2D is now due, 
presumably to the effects of the $ source-terms in the 
transport equation. Indeed, those terms are expected to 
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contribute as an effective "negative viscosity" , primarily 
due to the pressure-Hessian rotating small-scale strain 
matrices relative to the large-scale strain [HI]. Note that, 
strictly speaking, this is probably also true in 3D, so that 
the matrix-exponential closures are likely to be overly dis- 
sipative in every dimension. The main effect of the gra- 
dient stretching terms — which is a tendency to forward 
cascade, or positive eddy-viscosity — is well captured by 
the matrix-exponential closure in any dimension, but the 
additional, more subtle physics of the remaining terms in 
Eq. |4] are most likely not well-represented by the simple 
assumption of isotropy. 

As has been cautioned several times in this paper, the 
simplified matrix-exponential closure as written in Eq. 
1101 employs the drastic approximation of entirely omit- 
ting the pressure-strain correlation and other 'nonlin- 
ear scrambling' terms. But unlike eddy-viscosity based 
closure assumptions, this expression can be derived di- 
rectly from a relevant fluid dynamical equation, namely 
the stress transport equation (with only the production 
term), and using physically motivated and straightfor- 
ward assumptions about Lagrangian decorrelation and 
upstream isotropy. A similar result is obtained using an 
Eulcrian-Lagrangian change of variables when the stress 



is expressed in terms of subgrid-scalc velocity gradients. 
Perhaps it can be expected that casting this new light on 
the closure problem improves our understanding of this 
long-standing problem. 

Finally, we remark that many transport equations for 
turbulence moments have a basic structure similar to Eq. 
131 including two production terms involving the veloc- 
ity gradient and its transpose. Examples include higher- 
order moments of velocity, the spectral tensor encoun- 
tered in Rapid Distortion Theory calculations, etc... The 
formal solution in terms of matrix exponentials provides 
new possibilities of calculation and insights into the un- 
derlying physics. 
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